clear

global workingdata "working data"
global results "results"

//global doctab "male agenew1 pcertificate apcertificate rcertificate othertitle"
//global quality "r_qep e_qep corrdiag pcorrdiag cpcorrdiag corrmana pcorrmana cpmana medi nummedi_cond corrmedi_cond antibiotic_cond refer corrrefer_cond totalfee"
//global phy2 "male agenew1  apcertificate rcertificate othertitle"

***********************************************
*** Table 1: No Need to Replicate
***********************************************


***********************************************
*** Table 2: Summary of quality outcomes
***********************************************
use "$workingdata/append.dta", clear

global doctab "agenew1 male title"
global quality "r_qep corrdiag cpmana refer  medi corrmedi_cond corrmedi antibiotic_cond antibiotic"

*** To obtain mean and sd for all variables
eststo clear
	qui eststo: estpost tabstat $quality if (case == "A" & rural == 1), statistics(sum mean) columns(statistics)
	qui eststo: estpost tabstat $quality if (case == "A" & county == 1), statistics(sum mean) columns(statistics)
	qui eststo: estpost tabstat $quality if case == "A" & MVC == 1, statistics(sum mean) columns(statistics)
	qui eststo: estpost tabstat $quality if case == "A" & CHC == 1, statistics(sum mean) columns(statistics)
	qui eststo: estpost tabstat $quality if case == "A" & online == 1, statistics(sum mean) columns(statistics)
	
	qui eststo: estpost tabstat $quality if (case == "D" & rural == 1), statistics(sum mean) columns(statistics)
	qui eststo: estpost tabstat $quality if case == "D" & county == 1, statistics(sum mean) columns(statistics)
	qui eststo: estpost tabstat $quality if case == "D" & MVC == 1, statistics(sum mean) columns(statistics)
	qui eststo: estpost tabstat $quality if case == "D" & online == 1, statistics(sum mean) columns(statistics)
	
	qui eststo: estpost tabstat $quality if (case == "T" & rural == 1), statistics(sum mean) columns(statistics)
	qui eststo: estpost tabstat $quality if case == "T" & county == 1, statistics(sum mean) columns(statistics)
	qui eststo: estpost tabstat $quality if case == "T" & MVC == 1, statistics(sum mean) columns(statistics)
	qui eststo: estpost tabstat $quality if case == "T" & online == 1, statistics(sum mean) columns(statistics)
		
	esttab ///
		using  "$results/Table2.csv", ///
		main(sum %9.0fc) aux(mean %9.3fc) nostar unstack nolabel ///
		onecell ///
		nogap noeqlines nonote nolines obslast ///
		replace 

*** To obtain median and IQR for r_qep, saved in a separate file
eststo clear
	eststo: estpost tabstat r_qep if (case == "A" & rural == 1), statistics(p25 p50 p75 iqr) columns(statistics)
	eststo: estpost tabstat r_qep if (case == "A" & county == 1), statistics(p25 p50 p75 iqr) columns(statistics)
	eststo: estpost tabstat r_qep if case == "A" & MVC == 1, statistics(p25 p50 p75 iqr) columns(statistics)
	eststo: estpost tabstat r_qep if case == "A" & CHC == 1, statistics(p25 p50 p75 iqr) columns(statistics)
	eststo: estpost tabstat r_qep if case == "A" & online == 1, statistics(p25 p50 p75 iqr) columns(statistics)
	
	eststo: estpost tabstat r_qep if (case == "D" & rural == 1), statistics(p25 p50 p75 iqr) columns(statistics)
	eststo: estpost tabstat r_qep if case == "D" & county == 1, statistics(p25 p50 p75 iqr) columns(statistics)
	eststo: estpost tabstat r_qep if case == "D" & MVC == 1, statistics(p25 p50 p75 iqr) columns(statistics)
	eststo: estpost tabstat r_qep if case == "D" & online == 1, statistics(p25 p50 p75 iqr) columns(statistics)
	
	eststo: estpost tabstat r_qep if (case == "T" & rural == 1), statistics(p25 p50 p75 iqr) columns(statistics)
	eststo: estpost tabstat r_qep if case == "T" & county == 1, statistics(p25 p50 p75 iqr) columns(statistics)
	eststo: estpost tabstat r_qep if case == "T" & MVC == 1, statistics(p25 p50 p75 iqr) columns(statistics)
	eststo: estpost tabstat r_qep if case == "T" & online == 1, statistics(p25 p50 p75 iqr) columns(statistics)
	
	esttab ///
		using  "$results/Table2_medianiqr.csv", ///
		main(sum %9.0fc) aux(mean %9.3fc) nostar unstack nolabel ///
		onecell ///
		nogap noeqlines nonote nolines obslast ///
		replace 

*** p-values mentioned in the text
gen ruralvsmigrant = (level == "rural") if  (level == "rural"| level_sub == "MVC") & case == "A"
ttest cpmana, by(ruralvsmigrant)  level(95) 

gen CHCvsrural = (level_sub == "CHC") if (level_sub == "CHC" | level == "rural") & case == "A"
ttest cpmana, by(CHCvsrural)  level(95) 
ttest refer, by(CHCvsrural)  level(95) 		

*******************************************************************
*** Table 3: Association between quality outcomes and provider types
*******************************************************************
clear
use "$workingdata/append.dta", clear
gen lnfee = ln(1 + totalfee)
encode level_sub1, gen(level_sub1_id)
eststo clear

local i=0
foreach out of varlist ///
				 r_qep  ///
				 {
			local i=`i'+1
			eststo ols`i': reg `out' i.county  i.MVC  i.CHC i.online	i.D i.T ,cluster(stratification)
	}
local i=0
foreach out of varlist ///
				  corrdiag cpmana refer  medi corrmedi_cond corrmedi antibiotic_cond antibiotic ///
				 {
			local i=`i'+1
			logit `out' i.county  i.MVC  i.CHC i.online i.D i.T , cluster(stratification) level(95)	
			eststo logit`i': margins,dydx(*) post
	}
esttab ols1 logit1 logit2 logit3 logit4 logit5 logit6  logit7 logit8 ///
		using "$results/Table3.csv", scsv ///
			cells(b(fmt(3)) p(fmt(4) par(( ))) ci(fmt(3) par([ , ]))) ///
			ar2(2)  keep( 1.county 1.MVC 1.CHC 1.online 1.D 1.T ) ///
		  	order( 1.county 1.MVC 1.CHC 1.online 1.D 1.T ) ///
		    nostar r2 ///
			prehead(sep=;) ///
			scalar(mean sd) sfmt(%9.3fc) ///
			nogaps mtitles ///
			replace


/* CI: open scsv as text, insert sep=; in the top line. Enter. Then open scsv in excel. */

********************************************************************************
*** Table 4: Association between physician characteristics and quality outcomes.
********************************************************************************
clear
use "$workingdata/append.dta", clear
gen lnfee = ln(1 + totalfee)

eststo clear
local i=0
foreach out of varlist ///
				 r_qep  ///
				 {
			local i = `i' + 1
			eststo ols`i': reg `out' i.male i.agenew2 i.apcertificate i.rcertificate i.othertitle i.D i.T , cluster(stratification)
	}
local i=0
foreach out of varlist ///
				 corrdiag cpmana refer  medi corrmedi_cond corrmedi antibiotic_cond antibiotic ///
				 {
			local i=`i'+1
			logit `out' i.male i.agenew2 i.apcertificate i.rcertificate i.othertitle i.D i.T , cluster(stratification) level(95)	
			eststo logit`i': margins,dydx(*) post
	}
esttab ols1 logit1 logit2 logit3 logit4 logit5 logit6  logit7 logit8 ///
		using "$results/Table4.csv", scsv ///
			cells(b(fmt(2)) p(fmt(4) par(( ))) ci(fmt(2) par([ , ]))) ///
			ar2(2)  keep( 1.male 1.agenew2 1.apcertificate 1.rcertificate 1.othertitle 1.D 1.T ) ///
		  	order( 1.male 1.agenew2 1.apcertificate 1.rcertificate 1.othertitle 1.D 1.T ) ///
		    nostar r2 ///
			prehead(sep=;) ///
			scalar(mean sd) sfmt(%9.3fc) ///
			nogaps mtitles ///
			replace

********************************************************************************
*** Figure 1: 
********************************************************************************
use "$workingdata/append.dta", clear

global doctab "male agenew1 pcertificate apcertificate rcertificate othertitle"
global quality "r_qep e_qep corrdiag pcorrdiag cpcorrdiag corrmana pcorrmana cpmana medi nummedi_cond corrmedi_cond antibiotic_cond refer corrrefer_cond totalfee"

label variable male "Male"
label variable agenew1 "Age: <50"
label variable pcertificate "Practicing Certificate"
label variable apcertificate "Assistant Practicing Certificate"
label variable rcertificate "Rural Physician Certificate"
label variable othertitle "Others"

eststo clear
eststo est1: estpost ci $doctab
eststo est2: estpost ci $doctab if level == "rural"
eststo est3: estpost ci $doctab if level_sub == "County"
eststo est4: estpost ci $doctab if level_sub == "MVC"
eststo est5: estpost ci $doctab if level_sub == "CHC"
eststo est6: estpost ci $doctab if level == "online"

coefplot (est1, msymbol(o) mcolor(red) ciopt(lcolor(red))) /// 
		 (est2, msymbol(sh) mcolor(blue) ciopt(lcolor(blud))) ///
		 (est3, keep(male agenew*) msymbol(th) mcolor(orange) ciopt(lcolor(orange))) /// 
		 (est4, msymbol(dh) mcolor(purple) ciopt(lcolor(purple))) ///
		 (est5, msymbol(X) mcolor(green) ciopt(lcolor(green))) ///
		 (est6, msymbol(oh) mcolor(gray) ciopt(lcolor(gray))) ///
         ,ci((lb[1] ub[1])) ///
		 xlabel(0(0.25)1,format(%03.1f)) ///
		 legend(order(2 "Overall" 4 "Rural" 6 "County" 8 "MVC" 10 "CHC" 12 "Online") cols(3)) ///
		 graphregion(fcolor(white)) bgcolor(white)
graph export "$results/Figure1.pdf",replace

*** p-values mentioned in the text
gen onlinevsrural = (level == "online") if level == "online" | level == "rural" 
ttest agenew1, by(onlinevsrural)  level(95) 

gen onlinevscounty = (level == "online") if level == "online" | level_sub == "County" 
ttest agenew1, by(onlinevscounty)  level(95) 

gen onlinevsmigrant = (level == "online") if level == "online" | level_sub == "MVC" 
ttest agenew1, by(onlinevsmigrant)  level(95) 

gen CHCvsonline = (level_sub == "CHC") if level_sub == "CHC" | level == "online"
ttest agenew1, by(CHCvsonline)  level(95) 

gen CHCvsrural = (level_sub == "CHC") if level_sub == "CHC" | level == "rural"
ttest pcertificate, by(CHCvsrural)  level(95) 

********************************************************************************
*** Appendix Table A1: No Need to Replicate
********************************************************************************


********************************************************************************
*** Appendix Table A2: Quality Outcomes across three disease cases
********************************************************************************
use "$workingdata/append.dta", clear

global doctab "agenew1 male title"
global quality "r_qep corrdiag cpmana refer  medi corrmedi_cond corrmedi antibiotic_cond antibiotic"

*** To obtain mean and sd for all variables

eststo clear
	qui eststo: estpost tabstat $quality if (case == "A" & rural == 1), statistics(sum mean) columns(statistics)
	qui eststo: estpost tabstat $quality if (case == "A" & county == 1), statistics(sum mean) columns(statistics)
	qui eststo: estpost tabstat $quality if case == "A" & MVC == 1, statistics(sum mean) columns(statistics)
	qui eststo: estpost tabstat $quality if case == "A" & CHC == 1, statistics(sum mean) columns(statistics)
	qui eststo: estpost tabstat $quality if case == "A" & online == 1, statistics(sum mean) columns(statistics)
	
	qui eststo: estpost tabstat $quality if (case == "D" & rural == 1), statistics(sum mean) columns(statistics)
	qui eststo: estpost tabstat $quality if case == "D" & county == 1, statistics(sum mean) columns(statistics)
	qui eststo: estpost tabstat $quality if case == "D" & MVC == 1, statistics(sum mean) columns(statistics)
	qui eststo: estpost tabstat $quality if case == "D" & online == 1, statistics(sum mean) columns(statistics)
	
	qui eststo: estpost tabstat $quality if (case == "T" & rural == 1), statistics(sum mean) columns(statistics)
	qui eststo: estpost tabstat $quality if case == "T" & county == 1, statistics(sum mean) columns(statistics)
	qui eststo: estpost tabstat $quality if case == "T" & MVC == 1, statistics(sum mean) columns(statistics)
	qui eststo: estpost tabstat $quality if case == "T" & online == 1, statistics(sum mean) columns(statistics)
	
	esttab ///
		using  "$results/TableA2.csv", ///
		main(sum %9.0fc) aux(mean %9.3fc) nostar unstack nolabel ///
		onecell ///
		nogap noeqlines nonote nolines obslast ///
		replace 

*** To obtain median and IQR for r_qep, saved in a separate file
eststo clear
	eststo: estpost tabstat r_qep if (case == "A" & rural == 1), statistics(p25 p50 p75 iqr) columns(statistics)
	eststo: estpost tabstat r_qep if (case == "A" & county == 1), statistics(p25 p50 p75 iqr) columns(statistics)
	eststo: estpost tabstat r_qep if case == "A" & MVC == 1, statistics(p25 p50 p75 iqr) columns(statistics)
	eststo: estpost tabstat r_qep if case == "A" & CHC == 1, statistics(p25 p50 p75 iqr) columns(statistics)
	eststo: estpost tabstat r_qep if case == "A" & online == 1, statistics(p25 p50 p75 iqr) columns(statistics)
	
	eststo: estpost tabstat r_qep if (case == "D" & rural == 1), statistics(p25 p50 p75 iqr) columns(statistics)
	eststo: estpost tabstat r_qep if case == "D" & county == 1, statistics(p25 p50 p75 iqr) columns(statistics)
	eststo: estpost tabstat r_qep if case == "D" & MVC == 1, statistics(p25 p50 p75 iqr) columns(statistics)
	eststo: estpost tabstat r_qep if case == "D" & online == 1, statistics(p25 p50 p75 iqr) columns(statistics)
	
	eststo: estpost tabstat r_qep if (case == "T" & rural == 1), statistics(p25 p50 p75 iqr) columns(statistics)
	eststo: estpost tabstat r_qep if case == "T" & county == 1, statistics(p25 p50 p75 iqr) columns(statistics)
	eststo: estpost tabstat r_qep if case == "T" & MVC == 1, statistics(p25 p50 p75 iqr) columns(statistics)
	eststo: estpost tabstat r_qep if case == "T" & online == 1, statistics(p25 p50 p75 iqr) columns(statistics)
	
esttab ///
		using  "$results/TableA2_medianiqr.csv", ///
		main(sum %9.0fc) aux(mean %9.3fc) nostar unstack nolabel ///
		onecell ///
		nogap noeqlines nonote nolines obslast ///
		replace 


*** p-values
gen ruralvsmigrant = (level == "rural") if  (level == "rural"| level_sub == "MVC") & case == "A"
ttest cpmana, by(ruralvsmigrant)  level(95) 

gen CHCvsrural = (level_sub == "CHC") if (level_sub == "CHC" | level == "rural") & case == "A"
ttest cpmana, by(CHCvsrural)  level(95) 
ttest refer, by(CHCvsrural)  level(95) 		
		
********************************************************************************
*** Table A3: Physician Characteristics
********************************************************************************
clear

use "$workingdata/append.dta", clear

global doctab "male agenew1 agenew2 pcertificate apcertificate rcertificate othertitle"

eststo clear
eststo: estpost tabstat $doctab, statistics(sum mean) columns(statistics)
eststo: estpost tabstat $doctab if level == "rural", statistics(sum mean) columns(statistics)
eststo: estpost tabstat $doctab if level_sub == "County", statistics(sum mean) columns(statistics)
eststo: estpost tabstat $doctab if level_sub == "MVC", statistics(sum mean) columns(statistics)
eststo: estpost tabstat $doctab if level_sub == "CHC", statistics(sum mean) columns(statistics)
eststo: estpost tabstat $doctab if level == "online", statistics(sum mean) columns(statistics)


esttab ///
	using "${results}/TableA3.csv", ///
	main(sum %9.0fc) aux(mean %9.3fc) nostar unstack nolabel ///
	nogap noeqlines nonote nolines obslast ///
	replace 

********************************************************************************
*** Appendix Table A4, Figure A1, Figure A2: No Need to Replicate
********************************************************************************
